Do f{R) theories matter? 
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We consider a modified action functional with a non-minimum coupling between the scalar cur- 
vature and the matter Lagrangian, and study its consequences on stellar equilibrium. Particular 
attention is paid to the validity of the Newtonian regime, and on the boundary and exterior match- 
ing conditions, as well as on the redefinition of the metric components. Comparison with solar 
observables is achieved through numerical analysis, and constraints on the non-minimum coupling 
are discussed. 
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I. INTRODUCTION 

Modern cosmology faces two outstanding challenges, 
namely the existence and nature of dark energy and dark 
matter. Many theories have been put forward to address 
both issues: for dark matter, several candidates are avail- 
able, such as weak-interacting particles (WIMPs) aris- 
ing from extensions to the Standard Model {e.g. axions, 
neutralinos), etc.; for dark energy, "quintessence" mod- 
els consider the slow- roll of a scalar field [1, 2], amongst 
other candidates; others suggest that the averaging of 
inhomogeneities at a cosmological scale may yield an ef- 
fective scalar field, thus accounting for the dark energy 
component of the Universe [3] . A possible unification of 
both "dark" components has also been suggested, resort- 
ing to a scalar field model [4] or an exotic equation of 
state, as featured by the so-called modified Chaplygin 
gas [5]. 

A different approach assumes that no extra energy 
content is needed and that the fundamental laws and 
tenets of gravitation may be incomplete, perhaps just a 
low-energy approximation; as a consequence, modifica- 
tions of the Friedmann equation to include higher order 
terms in the energy density p (see e.g. [6] and references 
therein) have been proposed or, at a more fundamental 
level, changes to the action functional. A rather straight 
forward approach lies in replacing the linear scalar cur- 
vature term in the Einstein-Hilbcrt action by a function 
of the scalar curvature, f{R)', alternatively, one could re- 
sort to other scalar invariants of the theory [7] . This has 
led to some success in replicating the accelerated expan- 
sion of the Universe, while comparison with its evolu- 
tion throughout the different ages has yielded some con- 
straints and exclusions to the form of f{R); perhaps the 
most well-known proposal of this type is the Starobin- 
sky inflationary model, where a quadratic term in the 
curvature is added to the usual linear form (plus cosmo- 
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logical constant), f{R) = i? — A -I- aR^ [8]. Solar system 
tests could also bring further insight, mostly arising from 
the parameterized post-Newtonian (PPN) metric coeffi- 
cients derived from this extension of general relativity 
(GR). However, some disagreement exists in the commu- 
nity, with some arguing that no changes are predicted 
at a post-Newtonian level (see e.g. [9] and references 
therein) ; amongst other considerations, this mostly stems 
from an approach based either in the more usual metric 
affine connection (that is, where the affine connection is 
taken a priori as depending on the metric), or in the 
so-called Palatini approach [10] (where both the met- 
ric and the affine connection are taken as independent 
variables). As an example of a clear phenomenological 
consequence of this extension of GR, it has been shown 
that f{R) = /o^" theories yield a gravitational poten- 
tial which displays an increasing, repulsive contribution, 
added to the Newtonian term [11]. 

Notwithstanding the significant literature on these 
f{R) models, few steps have been taken to address an- 
other interesting possibility: not only that the curvature 
is non-trivial in the Einstein- Hilbert Lagrangian, but also 
that the coupling between matter and geometry is not 
minimum; indeed, these are only implicitly related in the 
action functional, since one expects that covariantly in- 
variant terms in Cm should be constructed by contrac- 
tion with the metric {e.g. the kinetic term of a real scalar 
field, g^'^ (j) ^^(j) . A non-minimum coupling would imply 
that geometric quantities (such as the scalar invariants) 
would explicitly show in the action; asides from theo- 
retical elegance, this could have deep phenomenological 
implications: indeed, in regions where the curvature is 
high (which, in GR, are related to regions of high en- 
ergy density or pressure), the implications of such the- 
ory could deviate considerably from those predicted by 
Einstein's theory [12]. Related proposals have been put 
forward previously to address the problem of the accel- 
erated expansion of the Universe [13] and the existence 
of a cosmological constant [14]. 

In this sense, the immediate question, posed in the ti- 
tle of this work, is: what are the implications for the be- 
haviour of matter under such conditions? Some work has 
been put forward concerning this issue, namely changes 
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to geodetic behaviour [12], the possibihty of modelhng 
dark matter [15] and the violation of the highly con- 
strained equivalence principle [16]. Perhaps the most 
important consequence of the mentioned studies is that 
energy may no longer be covariantly conserved, that is, 
V'^T^i/ ^ 0, where T^^ is the energy-momentum tensor 
of matter; this occurs because, due to the presence of 
extra terms in the equations of motion, the Bianchi iden- 
tities no longer imply in the covariant conservation of the 
energy-momentum tensor. 

This work addresses what we believe is the natu- 
ral proving ground for a non-minimally coupled gravity 
model: regions where the density may be high enough, 
to evidence some deviation from GR, although moderate 
enough so that effects are still perturbative - a star. The 
results rely upon and expand the methodology followed 
by the authors in previous studies [17, 18]. This paper 
is divided in the following sections: first, we present the 
model upon which the subsequent work is based; then, we 
develop the equations of motion, aiming at the modified 
Tolman-Oppenheimer-Volkoff (TOV) equation, with due 
care taken regarding the validity of the Newtonian regime 
and resulting modified hydrostatic equilibrium equation; 
afterwards, we insert the polytropic equation of state into 
the latter, and compute the necessary observables; a nu- 
merical session follows, where profiles and bounds are 
computed for the relevant quantities; finally, a discussion 
of our results is presented. 



II. THE MODEL 

Following the discussion of the previous section, one 
postulates the following action for the theory [12], 



T 



2 5{^Cm) 



(3) 



As stated before, the Bianchi identities, V^G^i, = 
imply the non- (covariant) conservation law 



XF2 

l + A/2 



{g^..Cm-T^.)y^'R , (4) 



and, as expected, in the GR limit A — > 0, one recovers 
the conservation law V^T^,^ = 0. 



A. Scope of application 

It is our stated purpose to arrive at the modified form 
of the TOV equation, the relativistic version of the hy- 
drostatic equilibrium condition. Prom Eq. (2), it is 
clear that a full treatment of the equations of motion 
is unattainable unless some specific form for fi{R) and 
/2(i?) is provided. Furthermore, the presence of both the 
pure curvature and non-minimum coupling, respectively, 
still constitutes a daunting analytical challenge; hence, 
and since one is mainly interested in the relevance of 
the efi'ects within a high curvature and pressure medium, 
where /2 should overwhelm the modification of the pure 
curvature term, /i — 2kR (neglecting the contribution of 
the cosmological constant) , it appears sensible to discard 
the latter; a treatment of the standard f{R) scenario with 
/2 = may be found in Rcf. [19]. 

Mathematically, the chosen approximation reads as 



S = 



-h{R) + [l + Xf2{R)]C, 



(1) 



where fi{R) (with i = 1,2) are arbitrary functions of the 
scalar curvature R, Cm is the Lagrangian density of mat- 
ter and g is the metric determinant. For convenience, the 
contribution of the non-minimum coupling of /2 is gauged 
through the coupling constant A (which has dimensions 
[A] — [,f2]~^)- The standard Einstein-Hilbert action is 
recovered by taking f2 — and /i — 2k{R — 2A), where 
K = d^/lQnG and A is the cosmological constant (from 
now on, one works in a unit system where c = 1). 

Variation with respect to the metric yields the 
modified Einstein equations of motion, here arranged as 



{Fi + 2AF2.^m) Ruf — -j^fiQiiv — 

{U^,u - 9f,uU) {Fi + 2XF2jCm) + [1 + A/2] T^. 



(2) 



where one defines = V^i^Viy for convenience, as well 

as Fi{R) = ,f'{R), and omitted the argument. The mat- 
ter energy-momentum tensor is, as usually, defined by 



^\fi-KR\ « |A/2£™| , 
\Fi-2k\ « 2\XF2£m\ , 

□) {F2Cm)\ 



(5) 



where the second and third inequalities indicate that this 
regime stems not just from the comparison between con- 
tributions to the action functional, but also those in- 
volved in the modified Einstein equations (2); also, notice 
that the perturbative condition A/2 (i?) <C 1 has not been 
enforced yet. The last inequality is satisfied if the follow- 
ing stronger conditions hold, 



dFi 



dr 
d^F, 



< 2 



d{F2Cm) 



dr"^ 



< 2 



dr 

dHF2Cm) 



(6) 



dr^ 



This said, one considers the simplest form /2 = R; this 
linear coupling could arise from a first order expansion 

of a more general f2{R) function, in the weak field en- 
viron of the Sun. Also, one takes Lm = a natural 
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choice for the Lagrangian density of an ideal fluid [20] - 
characterised by the standard energy-momentum tensor 



Tfu. = {(> + p)u^u,, - pg^i. 



(7) 



where p is the fluid's pressure, p its matter-energy density 
and u its 4- velocity vector, with = {uq, 0), = 1, so 

that Wo = 9oo^- Inserting the above expressions for /2 
and Cm, the inequahties (5) become 



l/il « 2\iK + \p)R\ 
\Fi\ ^ 2\K + Xp\ , 
dFi 



dr 



« 2\Xp'{r)\ 



dr^ 



(8) 
(9) 

«2|Ap"(r)| ,(10) 



where the prime denotes derivation with respect to the 
radial coordinate. Clearly, if /i = 2kR Fi = 2k, these 
are trivially satisfied. Also, notice that f2 = R implies 
that [A] = [f^^] = = M-2. 

One must now ascertain the form of In 
Ref. [2] it is shown that acceptable models with 
fia = 2k{R - aR^-"") {a > 0, [a] = [i?™] = Af^™ and 
< TO < 1 or /i6 = 2k{R + aR^ - A) (with aA < 1 and 
[a] = [i?^^] = M~^) are cosmologically viable; the latter 
may arise from the renomalizability of the theory near 
the Planck scale, with a ~ M'^ and M = 10^^ GeV [2]. 
For the perturbative regime to be valid, one must have, 
for the /i = /la case 



{l + ap)Ri,,- -R{g^, + aTi,,)= (13) 

^(□/ni/ - 9iJ.un)p + ^^/^^^ ' 

where one defines the parameter a = \/k = IGttGX, with 
dimension [a] = M~'^; accordingly, both ap and ap are 
dimensionless quantities. In the above form, the physi- 
cal meaning of the proposed model is more transparent: 
aside from a pressure-dependent term on the r.h.s., the 
most interesting modification occurs on the l.h.s.; firstly, 
the contribution of the Riemann tensor is modified by 
a factor 1 + ap; secondly, the scalar curvature is cou- 
pled not only to the metric 5^,^, but also to the energy- 
momentum tensor Tjuy. One thus gets a clear picture of 
the matter-geometry interaction, which occurs via scalar- 
tensor combinations. 

By taking the trace of the above equation, one obtains 



R = 



-3aKlO(p + T _ 3p - p + SaKliJp 
' K[2 + a{T- 2p)] ~ k[2 + a{p - 5p)] 



(14) 



having inserted T = = p — 3p. Interestingly enough, 
in the "strong" a — > 00 regime it yields the asymptotic re- 
sult R = SOp/ (p—^p)'- in this regime, a varying, low den- 
sity can give origin to extremely high curvatures, while 
an almost uniform, even if a high density fluid might yield 
a vanishingly small curvature. 



lAI » aK 



1 



|A| ^ (1 — m)maK 
|A| (1 — m)maK 



1 



dR 



1 



<PR 

dr 



— -(1 + ^)- 



(11) 



' dR\ 
. dr ) 



R 



while, for the /i = fu case (which can be derived from 
the above inequalities, setting to = — 1,q;^— a and 
ignoring the cosmological constant term). 



C. Perturbative regime 

Foreseeing a later comparison with solar observables, 
which one knows are well predicted by GR, it is now 

assumed that the effect of /2 yields only perturbative 
corrections: in the approximation IA/2I = |Ai?| <C 1, the 
scalar curvature is given by i? « (3p — p)/k, resulting in 
|(A/k)(3p — p)\ = |a(3p — p)| ^ 1; anticipating the New- 
tonian approximation p <^ p, this yields \a\p <C 1 and 
\a\p ^ 1. Inserting Eq. (14) in the Einstein equation, 
one obtains, after some algebraic manipulation. 



|A| » aK 
|A| > 2aK 
|A| » 2aK 



1 dR 
p' dr 



1 



d^R 



dr^ 



B. Equation of motion 



(12) 



k[2 + a{p - 3p)]R^, = (15) 
(3p - p)g^„ + 2(1 - ap)Tf_,„ + aK(4n^,^ - gf,^n)p , 

keeping only first order terms in a. This shall be the 
main tool for the derivation of the TOV equation. 



III. STELLAR EQUILIBRIUM 

A. Static, spherical symmetric scenario 



With the above choice for /i and /2, the equation of 
motion becomes 



Since one is dealing with an ideal, spherically symmet- 
ric system, in which temporal variations are assumed to 
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occur only at the cosmological scale H^^, and hence neg- 
ligible at an astrophysical time scale, one considers the 
Birkhoff metric (in its anisotropic form) , given by the line 
element 



(16) 



with d^l = r'^{dO'^ + sin^d dcjy^), so that = 
sine e(^+'")/2^ 

One can resort to the expression of the Riemann tensor 
and Eq. (15) to obtain the following intermediary step 
(the full algebraic derivation is shown in Appendix A) : 



= P- (17) 



K[2 + a{p-3p)] ^ + e" 

f (P + 3p) + ^ Ue-^Ooo + Se-'^Drr + 2^]P 



B. Newtonian limit 

Before continuing the derivation, it is opportune to ad- 
dress the issue of the validity of the Newtonian regime; 
clearly, establishing this limit will enable many fruitful 
simplifications in the calculations ahead. In the stan- 
dard derivation of the hydrostatic equilibrium equation, 
this arises from a set of simplifications imposed on the 
relativistic TOV equation; the later reads, 



, ^ G [p{r) +p(r)]K(r) +47rp(r)r3] 
r2 l-2Gme{r)/r 



(22) 



and the Newtonian approximation is valid if the following 
inequalities are satisfied 

r > 2Gme{r) , p{r) > p(r) , me(r) > 47rp(r)r^ , (23) 
yielding the non-relativistic hydrostatic equation of state 



One now defines the parameter me, here called the ef- 
fective mass, to distinguish it from its identification with 
the gravitational mass of the unperturbed GR scenario, 
derived from the Schwarzschild metric; it is given by the 
usual expression 



e-" = 1- 



2Gme 



(18) 



which yields 



2Gmi 



(19) 



where the prime denotes differentiation with respect to 
the radial coordinate. Inserting the above in Eq. (17) 
and solving for m^, one gets 



P'{r) 



-G 



p{r)me{r) 



(24) 



Since is is assumed that the coupling between mat- 
ter and geometry will only produce a perturbative effect, 
leading to the redefinition of mass through Eq. (21), it 
is clear that no changes occur regarding the validity of 
the Newtonian approximation. Hence, foreseeing the ap- 
plication of the following results to the Sun, where such 
regime is valid, one may simplify the intermediate calcu- 
lations and, when convenient, insert the inequalities (23) 
in order to simplify o-dependent terms. By the same to- 
ken, terms involving the coupling between a-dependent 
quantities and any covariant derivatives may be evalu- 
ated by taking their Newtonian counterparts, 



□ooP ^ 



-v'p' 



(25) 



2 + a{p- 3p) 



(20) 



4G 



2 + a{p- 3p) 



OrrP ■■ 

UeeP ■■ 



■p" + ^(t'p'^ 
■■ e~^rp' Ri rp' 



Notice that this does not imply that one is directly de- 
riving the Newtonian hydrostatic equilibrium equation, 
since leading order terms will not be approximated. 



In the perturbative regime, one may take only first-order 
terms in ap and ap, obtaining 



C. Tolman-Oppenheimer-Volkoff equation 



ar 
8G 



1^ , P 3p2 



(21) 



Irr "T j P > 



using AGk = l/A-K. The above expression clearly shows 
the perturbation to the purely gravitational mass, de- 
fined by m' = Airr'^p. 



By following a procedure similar to the one leading 
to Eq. (17) (depicted in Appendix A), one obtains the 
following equation: 



K[2 + a{p- 3p)] 



2Gme\ v' 2Gme 



(26) 



P + f (P - P) + f (se-'^Doo + 5e-n., - 2^\ p . 
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Substituting by the expressions for the covariant deriva- 
tives, one gets 



{27) 



5(1- 

Gmi 



2Gm, 



2GmA 3z/' 



1|2 

r 



Girif 



P 



One now introduces the approximations discussed in the 
previous subsection. Also, the approximation ~ 
m'g = Anpr^ is taken, since the perturbative corrections 
would produce second order terms in a. One obtains, 
after a little algebra. 



k{2 + ap) ^— 

/ a \ an 



2Gm, 



(28) 



bp" 



Gm'^ - 
5 — 



3i/ 
2 



Solving for u' , one gets 



V' 

2" 



^ (2 + gp) me + 4 (2 - gp) 7rpr3 
r{r-2Gme)[2 + a{p+lrp')] 



p" + {bnGpr - ^) p' 



5^// 
4 



r - 2Gme 2 + a (p + f rp') 
When g = 0, one recovers the standard expression 



u' (5 + '^T^pr^ 
"2 ~ r2 - 2Gmer 



(30) 



Linearizing with respect to a yields 



v' m.,, + 4:Trpr^ 

2- = ^ r^-2Gm.r +"^(^'^^ ' 



(31) 



where the function h{p, p) is defined through 

h{p,p)= (32) 
5,, /57r^9 3 Gm. 1 \ , 

after considering the inequality r ^ 2Gme and also 
birGpr"^ /2 <^ 1; taking the Sun's maximum, central, 
density pc = 1.622 x 10^ kg/m^, and the radius Rq — 
6.955 X 10« m, one gets b'KGpr'^/2c^ = 4.57 x IQ-^ < 1. 

Before continuing, one can rewrite the expression for 
the gravitational mass, obtained before, but imposing the 
limit e~'^ « 1, v' a' ~ 0, which yields 



□ooP~0 , UrrP^p" , BeeP^rp' . (33) 
This approximation, together with p <^ p, implies that 



m'e « 47rrV + ar^ + ^ - 2^ 

Similarly, one obtains 

-p-3aK{p" + ^p') 



(34) 



k{2 + ap) 



(35) 



2k 



3\p" + -p' 

T 



2k 



Although the expression for the scalar curvature R is 
somewhat involved, one can approximate its derivative by 
the unperturbed value R' = —p' /2k, losing only terms of 
order O(g^); indeed, one could write 



R' = 
\R' = 
with 



-^ + aF{p,p)^ 
-^+aXF{p,p) = 



(36) 



F{P,P) 



-2^+« . 



F{P,P) = - 2 



3 p" + 



(37) 



To write the modified equation of hydrostatic equilib- 
rium, one now resorts to the non- (covariant )conservation 
of the energy-momentum tensor, Eq. (4) . The expression 
for the covariant derivative is purely geometric; in order 
to derive an expression for p', one aims at the p = r com- 
ponent of the above equation; with our choice of /2 = i? 
and Cm = P, one gets the modified TOV equation 



/ z^'/ ^ 2A / 2GmA „, , , 
^ + 2^^ + ^^=iTxrV-^)^^^''^ 

-^{P + P)-P' -Y:^R' ^P' +app' , 

dropping higher order terms and using R' = —p' /2k. 

Thus, one finally obtains a set of three differential 
equations for the problem at hand: 



P - 



2XpR' = --{p + p) , 



— = G 



— 2Gmer 



3p" 




8G ' 


4Gr 






+ g 





A-nGpp I ~ 1" 



(39) 
(40) 

(41) 
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Substituting Eq. (39) into (41), after some algebra one 
gets the modified TOV equation, 



p' + G{p+p) 
'5 



p" - AnGpp 



me + Anpr 
— 2Gmer 
P' 



(42) 



- J J P + Pp' 



This yields the non-relativistic hydrostatic equilibrium 
equation, 



, , Gruep 
P + — 5— = a 



-p" - AwGpp 



r-~-^ ] P + Pp' 



(43) 



D. Polytropic equation of state 



Realistic stellar models rely on four differential equa- 
tions, together with appropriate definitions [21]; aside 
from the mass conservation condition (40) and the hy- 
drostatic equilibrium equation (42) or (43), these express 
energy conservation and transport, through 



and 



dL{r) 
dr 



dT{r) 
dr 



L{r) 



(44) 



(45) 



respectively. In the above, L{r) is the energy flow across 
a sphere of radius r, e is the energy generation rate per 
mass unit, and Ac is the conductivity coefficient. For 
given e and Ac, which account for the processes ongoing 
inside the star, one is left with five unknowns, m'^, p, p, 
T and L: an additional relation is needed, in the form of 
a suitable equation of state p = p{p). 

Many candidate equations of state and solar models are 
available, with varying degrees of sofistication, account- 
ing for effects such as chemical composition, solar mat- 
ter mixing, discontinuities between layers, heat difusion, 
etc.; two outstanding examples are the Mihalas-Hummer- 
Dappen and OPAL equations of state [22]. However, 
solving the above set of differential equations with a re- 
alistic equation of state requires heavy-duty numerical 
integration with complex code designs, and an analysis 
of the perturbations induced by a non- minimum coupling 
between matter and curvature is beyond the scope of the 
present study. 

Instead, one may resort to a very simplistic assump- 
tion, the so-called polytropic equation of state. This 
is commonly given hy p = Kp^^^^^^^, where K is the 
polytropic constant, po is the baryonic mass density and 
n is the polytropic index [21, 23, 24]. The polytropic 



index n interpolates between the basic thcrmodynami- 
cal processes: n = — 1 for isobaric, n = for isomet- 
ric, and an infinite polj^tropic index n for an isothermal 
sphere. Adiabatic processes yield n = 1/(7a — 1); with 
lA = Cp/cv the adiabatic coefficient. Several crude ap- 
proximations to relevant astrophysical systems are also 
obtained: n = 3/2 may model the degenerate star 
cores found in giant (gaseous) planets, white or brown 
dwarfs, and red giants; a polytropic index n = 5 yields a 
boundless system (that is, with non- vanishing density ev- 
erywhere), which was taken by Schuster as the first can- 
didate for a stellar system. Finally, a polytropic equation 
of state was used by Eddington in his proposal for the 
first solar model, with n = 3; naturally, it does not offer 
an accurate description of the solar interior, and has been 
deprecated by the following advancements. 

Nonetheless, the use of a polytropic equation of state 
is still of interest, due to its simplicity and ease of ma- 
nipulation, which render it a valuable tool in more the- 
oretically driven studies, as is the present case (as an 
example, the generalized Chaplygin gas is a polytrope 
with index n = — 1/(1 -|- a) [18]). Clearly, more realis- 
tic assumptions regarding the structure of the Sun would 
improve the final results; also, the procedure outlined in 
this work could also be applied to more exotic bodies, 
either through the use of an adequate polytropic index, 
or via a more realistic equation of state, possibly yielding 
a more stringent constraint on the coupling parameter A. 

Recall that the effective mass rUe is defined in terms 
of the energy density p, which appears in the energy- 
momentum tensor; these two quantities are related 
through p = pb + np, yielding the relation 



n/(n+l) 



+ np 



(46) 



However, since one is interested in probing the Newtonian 
regime which occurs in the Sun, the condition p <C p 
holds; therefore, p ~ ps, and one may take the form 
p = for the equation of state [25]. 

In order to transform the modified hydrostatic equilib- 
rium Eq. (43) into a differential equation with a single 
variable, one first rewrites Eq. (43) as 



-p" - AttGpp 



r- ] P + PP 



= -AnGp-a ( ^ + ^ - 27rGp2 



(47) 



and inserts the polytropic equation of state, p = 
written as p = PcO"{0 and p = Pc6'"+H?), 
with ^ — r/ro a dimcnsionlcss variable and r§ = 
{n + l)pc/4:nGp'^; as stated before, pc = 1.622 x 
10^ kg/m^ is the central density, and Pc = 2.48 x 10^^ Pa 
is the central pressure. One obtains the perturbed Lane- 
Emden equation for the function 9{^): 
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1 

e 



-^]-Nr, 



nn+l 



^ 3n- 1 



9' 4(n+l) 
I — 



1 + Aj- 



4^ 2 



i48l 



where the prime now denotes derivation with respect to 

the dimensionless radial coordinate ^, and one has de- 
fined Ac = apc and Nc = Pc/Pc — 1-7 x 10~^, for con- 
venience. Obviously, setting Ac = one recovers the 
unperturbed Lane-Emden equation [21] 



(49) 



One may try to solve Eq. (48) analytically around 
^ = 0, and compare with the solution of the unperturbed 
equation, given (in the vicinity of ^ = 0) by [21] 



6^ 120^ 

This calculation (outlined in Appendix B) yields 



(50) 



E. Perturbative solution 

Inspection of Eq. (48) shows that it is a third-order 
equation, while the unperturbed Lane-Emden equation is 

only second-order. Thus, the initial conditions 6{0) = 1 
and ^'(0) = alone are not sufficient. Indeed, an extra 
condition regarding the initial behaviour of the second 
derivative must be provided; this is derived from the se- 
ries expansion of 9 around ^ = taken before, that is 



In order to numerically solve the perturbed Lane- 
Emden equation, one rewrites Eq. (48) (after a long 
derivation) as 



with 



and 



A=Ul-Ac'';''\ 
6\ 12(n-hl) 



n / , 1139-|-59n 

B = 1 - Ac 

120 V 60 n+1 



(51) 



(52) 



(53) 



neglecting the Nc term and taking the limit |^c| 
1. Taking Ac = gives back the unperturbed values 
A — 1/6 and B = n/12G, as expected. 

For n = (a constant density model), one finds 



»({).i-l(i-f..U^ 



(54) 



For n = 3 (the first proposed model for the Sun), one 
obtains 



,(0. 1- 1 (l - 11^,) «^ + _L(i _»,.),«. (55) 

Finally, for n ^ oo (an isothermal sphere), one gets 



24 



5Ac^ 

nn+l 



9" + 



q/3 



Nc 



n{n — 1) 
7n 



1 9' 



(58) 



.3n + 2 9" 



n+1 C n+1 



+ 



4 



- + 2(2n + l)Nce' 



9' 



3^„^83n + 2.' 



5 n-l-1 ^ 



Clearly, when Ac — > 0, the first term blows up unless it is 
compensated by the condition 9" + 29' = —0", which 
is precisely the unperturbed Lane-Emden equation. By 
the same token, the first term expresses the deviation 
from the unperturbed case, and should be of order A^ 
therefore cancelling out the divergence. 

However, implementing the above third order differ- 
ential equation proves too computationally demanding; 
thus, one must first approach the perturbed Lane-Emden 
equation and, given that one is searching for a perturba- 
tion, expand it in terms of ^ = 9q{1 + Ac5), where 
is the solution to the unperturbed Lane-Emden equation 
and (5 is the (relative) perturbation, obeying |^c|<^ ^ 1- 
Thus, one may write, for n > 0, 



r = ^^(1 + Ac6T « + nAc5) 



(59) 



1 



23 



0(e)«i-- 1--A,U^ + - 1 



12 



1 



40 



1837 

Tso 



Ac e 



(56) 



If this expansion is introduced in the perturbed Lane- 
Emden equation, and considering only terms of order Ac, 
one obtains a differential equation for (5, with 6*0 as source: 



8 



S" + 2 + i) 5' + (n - l)9r'S = ^-^^0^-% (60) 
+(2n + l)N^^er-X + ^^C"' + 3A^c^o' 



:2n 



4(n+l) 

5n(n-l) 3 ,3 n(3n + 7) 



4(n + 1) 



2/1/2 

t^o 



1 orx 



having eliminated the second derivative of the unper- 
turbed solution through Eq. (49). From Eq. (51), one 
concludes that this differential equation is supplemented 
by the initial conditions S{0) = and S'{Q) = 0. 

Note that 5 does not depend on A^. one must only 
find the unique 6 (for each and n): one aims at a 
simultaneous variation of the model's parameter Ac and 
the polytropic index n in the vicinity of the standard 
solar value n = 3, enabling the plotting of an exclusion 
graph in the (Ac,n) plane. For n = 3, the differential 
equation for S simplifies to 



S" + 2\^-^ + jjS' + 2916 = SN^e^o + 2^o + (61) 



a'_iRQ) = -2G 

<T+ = -2 
u'_{Rq) = 2 



Me - MiR, 







Rq{Rq~2GM^) ' 
GM 

Rq{Rq - 2GM) ' 

GMe 



(63) 



Rq{Rq - 2GMe 



^(^p"(i?0)i?o-p'(i?o)) , 



1/'^ = -f7+ = 2 



GM 



i?0(i?0 - 2GM) 



where the + and — subscripts indicate inner or outer 
boundary condition. If one identifies M = Mg, it follows 
that = 0, which implies 

We)^y(M = o, (64) 
2 Rq 

(after resorting to Eq. (21), with p{Rq) = 0), and 



-p"(i?o)i?0-p'(i?0) = O . (65) 
Clearly, both conditions hold, since 



F. Mass budget and matching conditions 

One now computes the deviation between the effective 

mass TTig, defined in Eq. (21), and the gravitational mass, 
defined by m' = Anr^p; some algebra yields 



m^-rrig = — 



-AcN^e^ 



3n e'i 



+ — + 



^0 ;6^) 



which is explicitly written in terms of Ac = ape, the 
model's dimensionless parameter, and = Vcl Pc-, the 
ratio that measures the validity of the Newtonian regime. 
One has also to eliminate through Eq. (49); one con- 
cludes that not only is the above mass difference small 
(due to the pcrturbativc regime. \Ac\ <§; 1), but that it is 
further suppressed by the factor Nc <^ 1. 

Clearly, the model's parameters should be adjusted to 
the known observables: the star's radius Rq and mass 
M. The issue of this identification is rather delicate; 
indeed, M should not be identified with the gravitational 
mass mg(i?0), but with the total effective mass me(i?0). 
Asides from the physical interpretation of the full energy 
content of the star affecting geodetic motion around it, 
one can also resort to the matching conditions of the 
inner metric with the outer Schwarzschild metric; indeed, 
writing Me EE TOe(i?0), Mg = mg{RQ), = m'^RQ) 
and Mg = m'g{RQ), one has 



p'{Rq) « 0'(ei)r(a) = o , (66) 

p"{Rq) (X e"{c,)0{^,y^ + n0'{(^)'9{^,y^-^ = o , 

where signals the star's boundary, through 6{£,i) = 0. 



G. A solution for the divergence problem 

Given the identification M = Me = rne(i?0), one can 
now deduce the perturbation to the central density pc, 
which is a model dependent parameter. However, before 
using Eq. (48) to evaluate the perturbation S and extract 
relevant quantities, it should be noticed that, by inspec- 
tion, it is clear that S diverges when 0q approaches zero. 
Hence, one cannot extend the perturbational approach 
to the full range of the star, and should instead deal with 
the full differential equation for the perturbed 0. 

In order to circumvent this issue, recall that there is 
a pronounced deviation between the predictions of the 
polytropic model and the realistic Standard Solar Model, 
for r > Rr = 0.713 Rq; this refiects the crossing from the 
radiative zone, where the n = 3 polytrope is a good ap- 
proximation for the Sun, and the convection zone, where 
this approach fails. Hence, the irregular behaviour of S 
near may be safely disregarded, since the fundamen- 
tal equation for is not valid there. Instead, one shall 
consider only the range < r < Rr or, equivalently, 
< ^ < = 0.713 ^1. In doing so, one is of course 
neglecting the contribution of the latter for all relevant 
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quantities: however, although the density and pressure 
are still significant at that point, both the polytropic and 
the Standard Solar Model show that up to 99.1% of the 
Sun's mass is located within the radiative zone. 

Also, it can be numerically shown that the boundary 
condition 9{^) = 0o(O(l + A,SiO) = does not shift 
significantly from the unperturbed ^?o(Ci) = case, so 
that no problem arises from neglecting any changes to 
the scale factor ro (which, recall, relates the dimension- 
less coordinate ^ with the physical distance to the center 
r). Furthermore, the matching of the inner and outer 
derivatives of the metric should not be taken as a real- 
istic condition, but merely a consistency check for the 
developed model. 



m 
m 

0.4 

& 

-0.a 



n = 2.8 
■r} = 3: 




0.4 
r/R 



H. Model-dependent parameters 



FIG. 1: Profiles for the unperturbed solution 9o (top) and 
absolute perturbation A = S9o (bottom), for 2.8 < n < 3.2. 



One proceeds with the calculation of the total mass of 
the Sun or, more accurately, the mass of the radiative 
zone: for the unperturbed case, one has (see Appendix 
C for the derivation of the following results) 



M 



-^TT—^PcoO'or 



PcO 



M 



(67) 



where one defines = ^o(^r) and pco is the unperturbed 
central density (that is, obtained from M, R and the 
numerical results for the unperturbed solution Oq). 
In the perturbed case, one obtains 



interpret this integral as the volume averaging of the ex- 
pression in brackets within the integrand, with the den- 
sity as weighting function, that is, ^" oc p. 

From the polytropic equation of state, one gets p oc 
T"+^, which enables writing 



Tc 







n+1 



(70) 



n6 - 



3n^ 
8 Oo 



10_ 

2^ 



M 



i?3 

-4TT-^Pc9or 



(68) 



A, 



Or JO 



2nn 



Notice that there are two Ac-dependent contributions: 
one arising from 9q and its derivatives, since including 
the perturbation 5 would only amount to second-order 
terms O(A^), and other involving the integral of 5. 

Since the total mass M does not change (only its in- 
terpretation as a purely gravitational mass or a sum of 
gravitational plus "active" components), one gets 



PcO 



(69) 



'ir^Or JO 



n6 ■ 



3n 9'o^ 



2$ 8 



d^ 



Clearly, taking Ac = yields pc = Pco- Also, notice 
that is not perturbed: this would reflect a change in 
(since = 0.713 ^i), which has already been ruled out. 
Also, notice that the denominator affecting the integral 



is equal to ^"^9'^^ = /o*"(C^6'^)' = - 9^^: thus, one can 



IV. NUMERICAL RESULTS 

The above results allow us to construct an exclusion plot 
for the central temperature on the (n, A^) parameter 
space, by imposing the constraint 1 — Tc/Tco < 0.06 (6%), 
the uncertainty of the central temperature of the Sun. 
Solving numerically the differential equation for 5 for a 
varying n, and computing all relevant quantities for vary- 
ing Ac, and defining the absolute perturbation A = ^o^, 
so that 6 = ^o(l + Ac5) = ^o + Ac^, one obtains the 
results shown in the Figures 1 to 7. 

First, notice that, from Fig. 1, the absolute perturba- 
tion A is fairly insensitive to the value of n; clearly, given 
that the profile of 9q for varying n differs more sharply as 
^ ^ (as seen in Fig. 1), a approximately constant A 
translates into a relative perturbation 5 that also exhibits 
this behaviour, as can be seen from Fig. 2. Furthermore, 
notice the indication of the divergence of the relative per- 
turbation 5 in Fig. 2 (which does not depend on n, since 
^0^0 and A — > const, yields 5 ^ oo, as previously 
discussed); as stated before, this divergence is avoided 
by dealing only with the radiative region r < Rr. 

Also, as can be seen from Fig. 2, the relative perturba- 
tion 5 peaks at |^|max ~ 1-1 (for n = 2.9, and a smaller 
value of |(5|maa; ~ 0.8 for n = 3); hence, the perturbative 
condition |Ac|(5 <C 1 translates to | Ad <C 1, leading to the 
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FIG. 2: Profile of the relative perturbation 5 for 2.8 < n < 
3.2. 



FIG. 5: Relative deviation of the central temperature 
Tc/Tco — 1 as a function of Ac and n. 





FIG. 3; Fraction of gravitational to total mass as a function 
of Ac, for 2.8 < n < 3.2. 



FIG. 6: Contour plot for the relative deviation of the central 
temperature Tc/Tco— 1 as a function of Ac and n, with contour 
lines of step 0.1%. 



0.02 



I 




FIG. 4: Relative deviation of the central temperature 
Tc/Tco - 1 as a function of Ac, for 2.8 < n < 3.2. 
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FIG. 7: Mass profiles, for n = 3 and Ac = -1, 0, 1. 
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chosen (extreme) interval for simulation —l<Ac<l. 
Hence, according to Figs. 4, 5 and 6, no relative devia- 
tion of the central temperature above the experimentally 
determined level of 6% is attained. However, the values 
found, of the order of 1%, indicate that any future refine- 
ment of the experimental error of Tc could yield a direct 
bound on the parameter Ac. 

The above results show that Ac is presently uncon- 
strained by solar observables, aside from the perturbative 
condition \Ac\ <C 1 ^ |A| <C n/pc; taking k = c^/16ttG = 
2.41 X 10^2 kg.m/s^ and pc = 1.622 x 10^ kg/w?, this 
yields |A| <C 1.48 x 10^^ j or, in natural units, |A| ^ 
4.24 X 10^'^ eV~^. However, this study assumes that the 
effects arising from the non-trivial matter-geometry cou- 
pling supersede those from the non-linear curvature term 
in the modified Hilbert- Einstein action (1), as expressed 
by the inequalities (8), (11) and (12). Having numerically 
determined the perturbative solution 5, one can now re- 
evaluate these conditions. For this, one first replaces R 
by its unperturbed value (with the Newtonian approach 
p), i? « — /9/2k, since corrections would amount to 
second-order terms in A. 

Considering the scaling laws, p = /c»c^"'(C) ^'^d p = 
PcS^^^iS,), and the inequalities (11) and (12), some al- 
gebra (sec Appendix D) yields, for the /i = fia = 
2k{R - aR^-^) case. 



|A| » 



|A| » 





m+l 




a 




/ 2Nce 




' 2k \"+' 


n 1 





(71) 



n+l \pc 



2Nr0 



while, for the /i = fib = 2k{R + aR?) case, 



A » 
A > 



2Nce ' 

n a 



(72) 



n + lNcO ' 
A > — C„,_i(0 



n + lNcO 
where one defines the quantity 




FIG. 8: Profiles of Ca,™, for m, = — 1, 0, 1, superimposed on 
profile of 6o and 6'q, with n = 3. 



Varying m in its domain < m < 1 (for the /i = fia 
case) and to = — 1 (for the /i = fib case); its profile 
is depicted in Figure 8. One concludes that |C„,m| < 
1.5, for both the f\ = fia case (0 < m < 1) as well as 
the /i = fib case (to = —1). This, together with the 
constraint 9{^) > 0.1 (as can be seen from Fig. 1) and 
the quantity Nc = Pc/ipcC^) = 1-7 x 10~^ yields, for the 
fi = fia case, 

|A| > 2.8 X 10^3(g_48 X 1036)'"(l - m)ma eV-^^^+"'^ , 

(75) 

while, for the f\ = fu case, 



|A| > 6.6 X l{fa . 



(76) 



Using the previously discussed value, a = (10^^ GeV)^'^ 
, which arises from the Planck-scale renormalization of 
the theory, one gets 



|A| > 6.6 X lO"-*" eV-^ = (3.9 x 10** GeVy 



(77) 



Also, a recent paper has reported a relation between 
the coefficient of the quadratic term in /i and the PPN 
parameter 7 [26], 



(73) 



One can plot the function Cn^miO by evaluating by 
the unperturbed solution 9o, using Eq. (49) and taking 
n « 3. One obtains 



2k^ 



1-7 



27-1 



(78) 



Hence, the current experimental constraint 7— 1 = (2.1± 
2.3) X 10-^ [27] yields 



C'n,m(0 



(1 + 3m)0'o + pc 



0', (30'o + l0o)+0l 



a < 0.17Mp,2 = (3.0 X 10^^ GeV) 



(79) 



(74) where Mpi = ^/Kc/G = 1.2 x 10^'' GeV is the Planck 
mass. This yields 
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|A| > 7.4 X 10-51 eV-^ = (1.2 x 10^^ GeV) 



(80) 



would allow an upper bound to be placed on A; this close- 
ness between the prediction of the perturbative model 
and experiment also seems to validate the latter. 



Clearly, the /i = fu case does not impose any sig- 
nificant bound on A, with the above value laying many 
orders of magnitude below the upper bound, |A| <C 
4.24 X 10^^ eV~'^, arising from the perturbative treat- 
ment. 

Notice that the above inequalities (71), (72), (75), (76), 
(77) and (80) are not to be considered as restrictions on 
the parameters, but as conditions for the validity of the 
regime where the effects of the matter-geometry coupling 
supersede those of the the non-linear curvature term in 
Eq. (1). 



V. DISCUSSION 

In this work we have examined a model where the usual 
Einstein-Hilbert action is modified, not only by allowing 
non-linear curvature terms to appear, but by also en- 
abling an explicit, non- minimum coupling between the 
curvature and the Lagrangian of matter fields. We first 
assume that this matter-geometry coupling is linear in 
curvature and its effect allows the non-linear curvature 
terms to be neglected - an assertion that is qualified in 
the end of the study, for two relevant non-linear f{R) 
models. 

We then proceed and perform the necessary calcula- 
tions to ascertain the effect of the latter in the hydro- 
static equilibrium of an n ss 3 polytropc such as the Sun. 
We assume a perturbative regime to the usual Tolman- 
Oppenheimmer-Volkoff equation, and take into consider- 
ation the validity of the Newtonian regime in this mod- 
ified theory, as well as the redefinition of relevant quan- 
tities, which arc computed numerically and compared 
with Solar observables. This goal is achieved through the 
use of the (non-relativistic) polytropic equation of state 

p = ifp^'*'^^^": as stated before, this is a very simplistic 
description of the complex behaviour of solar matter, and 
has been superseded by much more elaborate equations 
of state; we use it in our study so to better illustrate 
the effects of the matter-curvature coupling on an eas- 
ily understandable, closed model. However, we should 
remark that this restrictive treatment, although advan- 
tageous from the theoretical point of view, might conceal 
some of the more intricate phenomenology found in stel- 
lar systems, which could affect our results. Clearly, a 
subsequent study should consider a more realistic solar 
model. 

The results allow us to conclude that no strong con- 
straints on the matter-geometry coupling from the com- 
parison between the model's predictions and current ex- 
perimental sensitivity, aside from the perturbative ap- 
proach Xf2{R) = XR ~ {X/k)Pc <^ 1, which yields 
|A| <C 4.24 X 10"^^ eV~'^. However, the numerically ob- 
tained results show that a slight increase in accuracy 



APPENDIX A 



Given the line element (16), one first writes 



-,00 p „rrr, _ „-a " 

J -Koo — g rirr — e 



(Al) 



Resorting to Eq. (15), one gets 



2k[2 + a{p - 3p)] (5"°i?oo - f^'Rrr) = (A2) 
2(1 - ap)(5°°Too - g^'Trr) + 4aK(5°°noo - ^"'□rr)^ ^ 



k[2 + a{p — 3p)]e' 



p + p - apip + p) + AaK{e '^□oo + e '^□rr)p , 
since the terms depending on the metric cancel out (that 

is, g°°goo - g'^'^grr = 0). 
Next, compute 



2K[2 + a{p-3p)]^ 



(A3) 



(3p - p)^ + 2(1 - ap)^ + ^{4006 - geeO)? ■ 



2k[2 + a{p - 3j3)] 



^ 1 
^ ^ • 2r 



p - 3p + 2(1 - ap)p + an -^Qe +n]p 



k[2 + a{p - 3p)] 



p-p- 2ap + OK -^Uee + □ P 



Adding Eqs. (A2) to (A3), one obtains Eq. (17) 



k[2 -t- a{p - 3p)] 



a' 1 



r 



(A4) 



- yCp + 3p) + Y ( 5e-"noo + Se-'^D^^ + ^0)9 j P 

with □ = e-'^Doo - er^Urr - '^UboIt'^. 

By subtracting Eqs. (A2) from (A3), and substituting 
Eq. (18), one gets Eq. (26) 



K [2 + a (p - 3p)] 



1 - 



2GtoA v' 2Gm 
r J r 



■ p (A5) 



+Y(p - P) + Y ( Se-'^Doo + 5e-"a^ - ^neej P 



13 



APPENDIX B one writes e{^) w l-A^'^+B^^. Inserting this in Eq. (48) 

yields, after some algebra and expanding up to fourth 
In order to obtain an approximation to the solution of order in ^, 
the Lane-Emden equation (48) in the vicinity of ^ = 0, 
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-1 + A. 



4(n+ 1) 
2 + 5A 



{A [ANcin + l)(2n + 1) + ^n(13 + 21n)] + 10B(13 + 21n)) 



(Bl) 











e+ 


An - Ac (^An 







Equating second and fourth-order terms, one gets 



and 



A^Ul-Ac^^, 
6\ 12(n+l) 



n / , 11 39 + 59n 

B « 1 - Ac 

120 I 60 n+1 



(B2) 



(B3) 



APPENDIX C 

In the unperturbed case, the total mass M is purely 
gravitational, and hence defined as usual: 

M = 4n r^p dr = 47rr^pco / f = (CI) 
Jo Jo 

-Anrlpco [^\ee'oy = -47rr3p,oe.'(^o)£=«. = 
Jo 

-An—^PcoOor . 

defining = ^o(^r); hence, the (unperturbed) central 
density pco, a model dependent parameter, is given by 



PcO 



(C2) 



In the perturbed case, one first uses Eq. (62) to define, 
for clarity. 



M 



r 

Jo 



i' dr 



(C4) 



47r dr + "^AcNc ^O^S dr = 

Jo ^ Jo 

Anrlpc ^ " eO"" d^ + ^^AcNcro ^ ' C'^o © d^ = 



1^ rir 



47r^Pe / e{Oo{l + AJ)y' d^ + 

,2 \ rir- 



Anrl^) A^Ncro [ ^6^0 d^ 

PcJ Jo 



Ri 



47r-f pe / eO^{l+nAc5) d^ + 



Anrlp^A, / e'^oQ dS, 



47r^p. (^-^eAr + Ac {n5 + 6) d^^ = 

-An^Pce'o, (i--^ t eo^ {nS + e) dA . 

\ 'ir^'Or Jo j 



Notice that there are two A^-dependent contributions: 

one arising from 6'o and its derivatives (embodied in G), 
and other involving the integral of b. 



e(0 = -^f 



10 

2e 



(C3) 



In the above, Qq is used instead of d, since 6 is coupled 
to Ac, so that including the perturbation b would only 
amount to second-order terms 0(^4^). Using the defini- 
tions Nc = Vcl Pc and Tq = (n -|- l)pc/47rGp^, one gets 



APPENDIX D 

By considering the scaling laws p = Pc&^iC) and p = 
Pc&^^^{0 and using Eq. (11) one gets, for the /i = f\a 

case, 
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|A|» 



|A| » 2™(1 - m)maK 



2k 

Pc 
1+m 



m+1 



P'ir) 



(Dl) 



pl+nipl 



n 



n+1 
|A|» 



2k\ 

Pc) 



m+l 



a(l — m)m 



2^(1 - m)maK}+"' 



2k\ 

Pc) 



m+l 



pl+ruplf 

a{l — m)m 
2A/'^6'"(i+™)+i 



p"{r) -{1 + m) 



p'Hr) 



n6l'2(0 + 6'6l"(0 



while, for the /i = /k, case, Eq. (12) yields 



where the prime denotes denotes differentiation with re- 
spect to ^ (the factors tq arising from changing derivation 
with respect to r to derivation with respect to ^ = r/ro 
cancel out). 
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